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We present a finite-element simulation tool for calculating light fields in 3D nano-optical devices. This allows to 
solve challenging problems on a standard personal computer. We present solutions to eigenvalue problems, like 
Bloch-type eigenvalues in photonic crystals and photonic crystal waveguides, and to scattering problems, like 
the transmission through finite photonic crystals. 

The discretization is based on unstructured tetrahedral grids with an adaptive grid refinement controlled 
and steered by an error-estimator. As ansatz functions we use higher order, vectorial elements (Nedelec, edge 
O ,. elements). For a fast convergence of the solution we make use of advanced multi-grid algorithms adapted for the 
O ' vectorial Maxwell's equations. 

Keywords: 3D photonic crystals, photonic crystal waveguides, finite-element method, Maxwell's equations, 
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-5 1- INTRODUCTION 

Photonic crystals (PhC's) are structures composed of different optical transparent materials with a spatially 
I ' periodic arrangement of the refractive index. 1,2 Propagating light with a wavelength of the order of the 
, periodicity length of the photonic crystal is significantly influenced by multiple interference effects. The most 
prominent effect is the opening of photonic bandgaps, in analogy to electronic bandgaps in semiconductor physics 
or atomic bandgaps in atom optics. Due to the fast progress in nano-fabrication technologies PhC's can be 
' manufactured with high accuracy. This allows for the miniaturization of optical components and a broad range 
of technological applications, like, e.g., in telecommunications. 3 The properties of light propagating in PhC's 
are in general critically dependent on system parameters. Therefore, the design of photonic crystal devices calls 
f^) , for simulation tools with high accuracy, speed and reliability. 

2. LIGHT PROPAGATION IN PHOTONIC CRYSTALS 

Light propagation in photonic crystals is governed by Maxwell's equations with the assumption of vanish- 
ed ' ing densities of free charges and currents. The dielectric coefficient e(x) and the permeability fi(x) are real, 
C^' positive and periodic, e (x) — e(x + a), n(x) = fi(x + a). Here a is any elementary vector of the crys- 
| tal lattice. 2 For given primitive lattice vectors ax, 02 and 03 the elementary cell O C K 3 is defined as 
fl = £ R 3 I x = aiai + a^o-^ + 0:303; < a\, ct2, 03 < 1 } • A time-harmonic ansatz with frequency u> and mag- 
netic field H(x, t) = e~ luJt H(x) leads to an eigenvalue equation for H(x) with the constraint that H(x) is 
■ divergence-free: 

Vx4rVxff(i)=wVf)iJ(i), V • fi(S )H (x) = 0, xen. (1) 

e(x) 

Similar equations are found for the electric field E(x,t) = e~ luJt E(x): 

V x — !— V x E(x) = uo 2 e(x)E(x), V • e(x )E(x) = 0, xgQ. (2) 
AW 
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The Bloch theorem applies for wave propagation in periodic media. Therefore we aim to find Bloch-type 
eigenmodes 2 to Equations (1), defined as 

H(x) = e lls u(x), u(x) = u(x + a). (3) 

where the Bloch wavevector k <E R 3 is chosen from the first Brillouin zone. A similar procedure yields the 
Bloch-type eigenmodes to Equations (2), however, in what follows we will concentrate on Equations (1). 

In order to reformulate Equations (1) and (3) we define the following functional spaces and sesquilinear forms: 

(a) The set of Bloch periodic smooth functions is defined as 

C2°(ft,C d ) = {io e C°° (fl,C d ) \ w (x + a) = e its w (£)} . 

The Sobolev space H^ (curl) is the closure of C2° (f2,C 3 ) with respect to the H (curl) -norm. The space Hi is 
defined accordingly. 

(b) The sesquilinear forms a : iJg (curl) x J2g (curl) — > C and b : H^ (curl) x H^ (curl) — > C are defined as 

a(w,v)= [ -(V x w) ■ (V x v) Ax, (4) 
Jo. £ 

b(w,v) = / [iw-vdx. (5) 



With this we get a weak formulation of Equations (1) and (3): 
Problem 1. Find uj 2 e R and H € H^ (curl) such that 

a(H,v) =uj 2 b(H,v) Vo £ ^ (curl) , (6) 

under the condition that 

b(H,Vp) = VpeHl. (7) 



The space H^ (curl) is the direct sum of the divergence-free subspace H~ (curl) and the subspace of gradient 
fields Vp, p G Hi. Hence, h e H^ (curl) can be decomposed as 

h = h^ + Vp 

(Hclmholtz decomposition), where p solves the equation 

J Vp- Vvdx = J h- Vvdx V v e Hi. 

3. FINITE ELEMENT DISCRETIZATION 

It is crucial to inherit the properties of the Helmholtz decomposition to the sub-spaces on the discrete level. 
Otherwise the discrete spectrum is polluted by many unphysical fields - called spurious modes - stemming from 
the space of gradient fields. 4 Using Nedelec's edge elements to discretize the space H% (curl) and standard 
Lagrange elements of the same order to discretize the space Hi gives a discrete counterpart to the Helmholtz 
decomposition and to the divergence condition. 5 

We denote the discrete subspaces as follows: W h £ C i?g (curl), V h £ C Hi. Bloch periodicity is enforced by a 
multiplication of basis functions associated with one of two corresponding periodic boundaries of the unit cell by 
the Bloch factor exp(ik ■ Si) (see Equation (3)). All interior basis functions remain unchanged. An alternative 
approach is discussed by Dobson ct al 6 : In this approach a wave equation is formulated for the periodic part of 
the magnetic field, u(x) (see Equation (3)). Modified finite element ansatz functions are constructed from the 

Sobolev space H ( curl + ik x ) . 



The discretized problem corresponding to Problem 1 reads as follows: 
Problem 2. Find uj 2 e K and H eW h ^ such that 

a(H,<j>)=u; 2 b(H,cl>) ¥^6^ (8) 

under the condition that 

b(H,Vp) = V P €V h j. (9) 

The finite element basis functions for W h g are denoted by d>j>, 1 < f < N c the basis functions for V h g by 
Vfe'i 1 < ^' < ATp- We expanding in if = ^Uj</>j. Inserting into Equation (8) yields the algebraic 

eigenvalue problem and the algebraic divergence condition 

Au = ABu (10) 
G h Bu = (11) 

with Ajj := afoifrfij), Bi.j = 6(0^0^) and dj defined by Vy>j — J2jGj,i4>j- The matrix A is hermitian, 
positive semidefinite and B is hermitian, positive definite. In the algebraic form the Helmholtz decomposition 
reads as u = u 1 - + Gp, where p solves the algebraic problem 

G h BG P = G fc Bu (12) 
Due to the locality of the finite element basis functions all matrices are sparse. 

4. NUMERICAL SOLUTION OF THE EIGENVALUE EQUATION 

To solve the algebraic equation (10) we use a preconditioned Dohler's method 7 which is based on minimizing 
the Raylcigh quotient. To avoid that the iteration tumbles into the non-physical kernel of the (curl)-operator we 
project the iterates onto the divergence-free subspace. 

We use multi-level algorithms 8 for preconditioning as well for performing the Helmholtz decomposition (12). 
This is similar to the implementation by Hiptmair et al. 9 

With this, the computational time and the memory requirements grow linearly with the number of un- 
knowns. 10 Furthermore, we have implemented a residuum-based error estimator 11 and adaptive mesh refinement 
for the precise determination of localized modes (see chapter 6). As FE ansatz functions, we typically choose 
edge elements of quadratic order. 4 

5. BAND STRUCTURES OF 3D PHOTONIC CRYSTALS 

A model problem for 3D photonic crystals are so-called scaffold structures. 12 The geometry of a unit cell 
(sidelength o) of a simple cubic lattice is shown in Fig. 1. It consists of bars (width d — 0.25 a) of a transparent 
material with relative permittivity e r = 13 and a background with e r = 1 (e = e r eo, Eq: free space permittivity). 
For the calculation of the band structure the Bloch wavevector k is varied along symmetry lines of the Brillouin 
zone (cf. Dobson 12 ). 

The band structure for light propagating in the scaffold structure is shown in Fig. 2a. It exhibits a complete 
bandgap around the reduced frequency of Co = ua/{2-nc) ~ 0.4 which is indicated by the dotted horizontal 
lines in Fig. 2a. Table 1 shows the four lowest eigenvalues at the X-point(fc = (n/a, 0,0)) calculated on grids 
generated in 0, 1, rcsp. 2, uniform refinement steps from a coarse grid. In each uniform refinement step, every 
tetrahedron is subdivided into eight new tetrahedra. Shown are also the numbers of unknowns in the problem 
(number of ansatz functions in the finite element discretization) and typical computation times on a PC (Intel 
Pentium IV, 2.5 GHz). It can be seen that the computational effort grows linearly with the number of unknowns. 
Figure 2 (b) shows the dependence of the relative error of the four lowest eigenvalues flu^jv — ^i,q\/^i, q ) on the 
number of ansatz functions in the expansion of the eigenfunctions (number of unknowns). Here, w^jv is the i th 
eigenfrequency of the discrete solution with N unknowns, is the i th eigenfrequency of the quasi-exact solution 
obtained from a calculation on a finite-element grid with N = 1764048 unknowns. 



Figure 1. Unit cell of a 3D photonic crystal {scaffold). Bars with quadratic cross-sections intersect and form a 3D 
structure, periodic boundary conditions apply to all pairs of opposing faces. 




6. CALCULATION OF DEFECT MODES USING ADAPTIVE GRID REFINEMENT 

Light at a frequency inside the bandgap of a photonic crystal can be "trapped" inside defects of the structure. 1 
This enables the construction of, e.g., waveguides (line defects) and micro-cavities (point defects). 

Figure 3 (a) shows the geometry of a 2D photonic crystal with a point defect (a missing pore in the center). 
It consists of a hexagonal lattice of air holes with a radius of r = OA a in a material with a relative electric 
permittivity of e r = 13. A corresponding coarse triangular FE grid is shown in Figure 3 (b). Please note that 
circular air pores are approximated by polygons. In the coarse grid shown in Figure 3 (a) the pores close to 
the center are approximated to a higher accuracy than pores in the outer regions. Obviously, when refining the 
coarse grid in order to get a discrete solution which is closer to the solution of Problem 1 in some norm, the best 
strategy will not be to refine the grid uniformly, but to refine it in certain regions onlz. For this adaptive grid 
refinement we have implemented a residuum-based error estimator. 11 

Figure 4 (a) shows the modulus of the magnetic field for the lowest-frequency trapped eigenmode, computed 
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Table 1. First eigenvalues of eigenmodes of the scaffold structure at k = X. Shown are the step number, the number of 
degrees of freedom of the problem, the CPU time (run on a standard PC), and the reduced frequencies of the four lowest 
eigenmodes. 




(a) (b) 
Figure 3. Geometry (a) and coarse FE mesh (b) of a 2D photonic crystal structure with a central point defect. 

with adaptive refinement of the FE mesh. Figure 4 (b) shows the convergence of the eigenvalue corresponding to 
this solution towards the eigenvalue of a quasi-exact solution for adaptive grid refinement and for uniform grid 
refinement. Obviously, adaptive grid refinement is especially useful when the sought solutions are geometrically 
localized, or when the geometry exhibits sharp features, like discontinuities in the refractive index distribution. 
In this example, the use of the error estimator and adaptive refinement yields an order of magnitude in the 
accuracy of the error for a number of unknowns of N ~ 10 5 . 

7. PHOTONIC CRYSTAL SLAB WAVEGUIDE 

Photonic crystal waveguides are promising candidates for a range of applications of micro- or nano-optical 
elements like dispersion compensators or input lines for further miniaturized elements. We examine a waveguiding 
structure composed of a slab waveguide (confinement of the light in z-direction by a high-index guiding layer) 
combined with a 2D hexagonal array of air holes (PhC). 13 The waveguide is formed by a defect row of missing air 
holes in T-K-direction. Therefore, in certain wavelength ranges, the light is confined vertically by total internal 
reflection in the guiding layer and horizontally by the photonic bandgap due to the 2D photonic crystal. We 
consider a guiding layer of height z = 200nm and refractive index n = 3.4 with a substrate and superstrate of 
z = 900nm each and refractive index n = 1.45, and six rows of pores with refractive index n = 1.0 on each side 
of the waveguide. These parameters correspond to the material system of air pores in a Si02 — Si — S1O2 slab 
structure. The pore radius is r ~ 0.36a, the lattice vectors have a length of a — 532 nm. 

Due to the symmetry of the problem it is sufficient to restrict the computational domain to one quarter of a 
unit cell. This reduced computational domain is shown in Figure 5(a). Here, the guiding layer is colored in dark 
gray, the substrate in gray, and the air pores in light gray. Mirror symmetries are applied to the upper plane (z: 
center of the guiding layer) and on the left (x = 0). Periodic boundary conditions are applied to the front and 
back planes. At the left of the computational domain a Wl waveguide (i.e., width W = 1.0 a) is formed by a 
missing row of air pores. 

Fig. 5(b) shows some of the tctrahcdral elements in the discretization of the geometry. A typical coarse grid 
for this problem consists of about 10 4 tetrahedra. 
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Figure 4. (a) Distribution of the magnetic field intensity (\H(x, y)\) for the lowest-frequency bound state at the point 
defect, (b) Comparison of the convergence of the eigenfrequency of the lowest frequency bound state towards a quasi-exact 
solution for adaptive and uniform refinement of the FE mesh. 

Solving equation (10) for this problem on a PC (Intel Pentium IV, 2.5 MHz, 2Gbyte RAM) typically takes 
our FEM code a time of about 2 min and delivers the eigenvalue and the complex vector field for given Bloch 
wavevector k. Figure 6 shows a specific guided mode in this structure. Part (a) of this figure shows the 
amplitude of the magnetic field in a gray scale representation. The plotted solution has been calculated for a 
Bloch wavevector of ka/(2n) = 0.23 and corresponds to an eigenvalue of u}a/(2irc) ~ 0.326, which lies inside 
the first bandgap of the (2D) photonic crystal. As can be seen from the figure, the light field is localized in the 
high index guiding layer and in the region of the missing pore. White lines indicate equally spaced iso-intensity 
surfaces. Part (b) of the figure shows a cross-section through the same (vectorial) solution in the upper mirror 
plane (z — const.). In this plane the electric field vectors are oriented in the x — y-plane, this solution corresponds 
to a TE-like mode. 

8. TRANSMISSION THROUGH A FINITE PHOTONIC CRYSTAL 

In order to simulate the transmission of an incident light field Ui„ with a frequency to through a photonic 
crystal of finite width we perform scattering calculations. In this case we have to seek a solution u which fulfills 
Equation (1) resp. (2) for the given frequency lo on the computational domain Q with the following boundary 
condition: (a) the field on the boundary can be written as a superposition of incident and scattered light field: 





Figure 5. (a) Geometry of the reduced unit cell of a Wl waveguide (dark gray: guiding layer, gray: substrate, light gray: 
cylindrical air pores), (b) Visualization of the tetrahedral discretization of the geometry. 
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Figure 6. (a) Magnitude of the electric field of a guided TE like mode in a gray scale representation (black: high 
intensity) . Iso-Intensity surfaces are indicated with white lines, (b) Cross section of the E field distribution in the upper 
mirror plane (z = const.). Additionally, the material distribution is indicated. 



u = Ui n + u sc , (b) the scattered field is purely outgoing {radiation condition), which gives a condition for the 
outward normal derivative of u, d v u sc . We obtain a relation between d„u sc and u sc by using a modified type of 
perfectly matched layer boundary conditions, which also allows to treat certain types of inhomogeneous exterior 
domains. 14 

In the following example we examine the light transmission through a finite 2D photonic crystal consisting 
of four rows of air pores in a high index material. Figure 7 shows the geometry of the problem: Air pores 
(n = 1.0) of radius r — 0.366 a in a high-index material (n = 3.4) form a finite hexagonal pattern. Periodic 
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Figure 7. Geometry of a 2D finite photonic crystal: Bloch-periodic boundary conditions apply to the boundaries (BC 
pi) and (BC p2), transparent boundary conditions apply to the boundaries (BC tl) and (BC t2), a plane wave is incident 
onto the boundary (BC tl). 



boundary conditions apply to two boundaries of the computational domain. Transparent boundary conditions to 
the exterior (which is assumed to be also filled by the high index material) are applied to the two other facets of 
the computational window. A plane wave with a freely chosen angle of incidence a is incident onto one boundary. 
In order to measure the field transmission through the domain, the field on the opposing facet (BC t2) detected. 

Figure 8 shows calculated field distributions for different frequencies Co = u>a/2irc = a/A of the incident plane 
wave (a = lOdeg). Each of these solutions has been calculated using adaptive mesh refinement, finite elements 
of quadratic order, and typically N ~ 4 • 10 4 unknowns. The total calculation time on a standard laptop (Intel 
Pentium IV, 2.0 GHz) for each solution amounts about 10 sec. The fields in Figure 8 (a)-(j) correspond to 
decreasing wavelength. It can easily seen how dramatically the transmission is changes with the wavelength: 
Certain wavelengths (b, g, h, j) correspond to (partial or full) band-gaps of the photonic crystal, where light 
transmission is suppressed. At other wavelengths the light is transmitted; resonance behavior / slow group 
velocities can be discovered by the observed increased field amplitudes in the regions between the air holes, see 
e.g. (f). 

We decompose the field at the facets of the computational window into Fourier series: 
u (y) — S^l-oo A n exp(i2irny/L). The relative transmission through the photonic crystal is then given by 

_ Z)|fc n |<|fc in | sm(k ni n)A n 
A 2 

where n is the normal vector on the end facet. Figure 9 (a) shows the transmission in dependence on the 
frequency of the incident light for an angle of incidence a = deg. Detailed structures corresponding to bands 
and band gaps can be observed. Figure 9 (b) shows the relative error of the lowest order Fourier coefficient of 
the transmitted light field in dependence on the number of ansatz functions. For TE light fields (\E z (x,y)\), 
the errors are lower due to the smoothness of the electric field in this case. However, even for TM light fields 
(\H z (x,y)\) very accurate transmission coefficients with errors in the 10~ 3 -range can be gained with rather low 
numbers of unknowns and in short computation times on standard PC's. 



(13) 




Figure 8. Light fields propagating through a finite 2D photonic crystal. Plotted is the magnitude of the magnetic field, 
\H z (x,y)\, in a gray scale representation. The direction of the plane wave is indicated by the vector which is incident 
under an angle of a = 10 deg (compare Fig. 7). The different plots correspond to incident plane waves with different 
frequencies: (a): d> = uja/2nc = 0.182, (b): u> = 0.303, (c): u> = 0.385, (d): Co = 0.4, (e): Cj = 0.5, (f): ui = 0.625, (g): 
uj = 0.667, (h): u> = 0.714, (i): ui = 0.741, (j): ui = 0.769. 

9. CONCLUSION 

In this paper we have presented an adaptive finite-element method solver for the computation of electromagnetic 
eigenmodes and scattered light fields. The convergence analysis of solutions for model problems shows the 
efficiency of the methods. Our solver has been shown to give very accurate solutions for typical problems arising 
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Figure 9. (a): Transmission of TE and TM light fields through a finite 2D photonic crystal for an angle of incidence 
of a — Odeg. (b): Convergence of the relative error of the zero order transmission for TE and TM light fields (for a 
frequency of u ~ 0.37 and uniform grid refinement). Cpu times for computations on a standard laptop are indicated 
(Intel Pentium IV, 2.0 GHz). 

in nano- and micro-optics - even challenging 3D design tasks can be tackled on standard personal computers. 
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